R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

#library(spdep)
#library(maptools)
#library(rgeos)
#library(GISTools)
library(sf)
## Warning: package 'sf' was built under R version 4.2.2
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
#library(sp)
#library(tmap)
#library(spatialreg)
#library(rgdal)
library(tidyverse)
## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(biscale)
## Warning: package 'biscale' was built under R version 4.2.2
#library(cowplot)
library(mapview)
## Warning: package 'mapview' was built under R version 4.2.2
setwd("C://Users//Kaifs//OneDrive//Documents//DHS individual//shpfile")

BJ <- st_read(dsn="BJGE71FL/BJGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `BJGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\BJGE71FL\BJGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 555 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 5.684342e-14 ymin: 5.684342e-14 xmax: 3.611158 ymax: 12.34528
## Geodetic CRS:  WGS 84
BU <- st_read(dsn="BUGE71FL/BUGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `BUGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\BUGE71FL\BUGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 554 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: -4.432215 xmax: 30.784 ymax: 0
## Geodetic CRS:  WGS 84
CM <- st_read(dsn="CMGE71FL/CMGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `CMGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\CMGE71FL\CMGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 430 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 9.200917 ymin: 2.028929 xmax: 16.06679 ymax: 12.76956
## Geodetic CRS:  WGS 84
GM <- st_read(dsn="GMGE81FL/GMGE81FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `GMGE81FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\GMGE81FL\GMGE81FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 280 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -16.80157 ymin: 13.14906 xmax: -13.84872 ymax: 13.80121
## Geodetic CRS:  WGS 84
GN <- st_read(dsn="GNGE71FL/GNGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `GNGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\GNGE71FL\GNGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 401 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -14.79391 ymin: 7.408167 xmax: -7.856225 ymax: 12.51171
## Geodetic CRS:  WGS 84
NG <- st_read(dsn="NGGE7BFL/NGGE7BFL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `NGGE7BFL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\NGGE7BFL\NGGE7BFL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1389 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 0 xmax: 13.66089 ymax: 13.76644
## Geodetic CRS:  WGS 84
LB <- st_read(dsn="LBGE7AFL/LBGE7AFL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `LBGE7AFL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\LBGE7AFL\LBGE7AFL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 325 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -11.41053 ymin: 0 xmax: 0 ymax: 8.419803
## Geodetic CRS:  WGS 84
ML <- st_read(dsn="MLGE7AFL/MLGE7AFL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `MLGE7AFL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\MLGE7AFL\MLGE7AFL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 345 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -11.92693 ymin: 5.684342e-14 xmax: 1.42239 ymax: 18.45378
## Geodetic CRS:  WGS 84
MR <- st_read(dsn="MRGE71FL/MRGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `MRGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\MRGE71FL\MRGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1200 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -17.06289 ymin: 0 xmax: 0 ymax: 25.42278
## Geodetic CRS:  WGS 84
RW <- st_read(dsn="RWGE81FL/RWGE81FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `RWGE81FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\RWGE81FL\RWGE81FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 500 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 28.89335 ymin: -2.8127 xmax: 30.82515 ymax: -1.078855
## Geodetic CRS:  WGS 84
SL <- st_read(dsn="SLGE7AFL/SLGE7AFL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `SLGE7AFL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\SLGE7AFL\SLGE7AFL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 576 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -13.28709 ymin: 0 xmax: 0 ymax: 9.976958
## Geodetic CRS:  WGS 84
ZM <- st_read(dsn="ZMGE71FL/ZMGE71FL.shp") %>% 
  dplyr::select("DHSCC","DHSYEAR","DHSCLUST","LATNUM","LONGNUM","geometry")
## Reading layer `ZMGE71FL' from data source 
##   `C:\Users\Kaifs\OneDrive\Documents\DHS individual\shpfile\ZMGE71FL\ZMGE71FL.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 545 features and 20 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -0.0048541 ymin: -17.94572 xmax: 33.60199 ymax: 0.00149906
## Geodetic CRS:  WGS 84
map <- rbind(BJ,BU,CM,GM,GN,NG,LB,ML,MR,RW,SL,ZM)


#label var cum_pre_hdevi   "cumulative abnormally high precipitation (m)"
#label var cum_pre_hdevi_E "cumulative abnormally high precipitation: early childhood(m)"
#label var cum_pre_hdevi_L "cumulative abnormally high precipitation: mid/late childhood (m)"

#label var cum_pre_ldevi   "cumulative abnormally low precipitation (m)"#
#label var cum_pre_ldevi_E "cumulative abnormally low precipitation: early childhood(m)"
#label var cum_pre_ldevi_L "cumulative abnormally low precipitation: mid/late (m)"

#//temperature
#label var cum_mtemp_hdevi   "cumulative abnormally high temperature (m)"
#label var cum_mtemp_hdevi_E "cumulative abnormally high temperature: early childhood(m)"
#label var cum_mtemp_hdevi_L "cumulative abnormally high temperature: mid/late childhood(m)"

#label var cum_mtemp_ldevi   "cumulative abnormally low precipitation (m)"
#label var cum_mtemp_ldevi_E "cumulative abnormally low precipitation: early childhood (m)"
#label var cum_mtemp_ldevi_L "cumulative abnormally low precipitation: mid/late childhood (m)"

env_africa<-read.csv("TenAfricacountries.csv")
env_africa$DHSCC <- env_africa$dhscc
env_africa$DHSCLUST <- env_africa$dhsclust

env_africa<-env_africa %>% 
  group_by(dhscc,dhsclust,dhsyear,DHSCLUST,DHSCC,all_population_count_2015,global_human_footprint,nightlights_composite) %>% 
  summarise(cum_pre_hdevi_m = mean(cum_pre_hdevi),
            cum_pre_hdevi_E_m = mean(cum_pre_hdevi_E),
            cum_pre_hdevi_L_m = mean(cum_pre_hdevi_L),
            cum_pre_ldevi_m = mean(cum_pre_ldevi),
            cum_pre_ldevi_E_m = mean(cum_pre_ldevi_E),
            cum_pre_ldevi_L_m = mean(cum_pre_ldevi_L),
            cum_mtemp_hdevi_m = mean(cum_mtemp_hdevi),
            cum_mtemp_hdevi_E_m = mean(cum_mtemp_hdevi_E),
            cum_mtemp_hdevi_L_m = mean(cum_mtemp_hdevi_L),
            cum_mtemp_ldevi_m = mean(cum_mtemp_ldevi),
            cum_mtemp_ldevi_E_m = mean(cum_mtemp_ldevi_E),
            cum_mtemp_ldevi_L_m = mean(cum_mtemp_ldevi_L),
            average_prim = mean(pri_com1)) %>%
  left_join(map, by=c("DHSCLUST","DHSCC"))
## `summarise()` has grouped output by 'dhscc', 'dhsclust', 'dhsyear', 'DHSCLUST',
## 'DHSCC', 'all_population_count_2015', 'global_human_footprint'. You can
## override using the `.groups` argument.
# right now try with point in time variable
#BJ_map <- unique( env_BJ[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
#  left_join(BJ, by="DHSCLUST")
#CM_map <- unique( env_CM[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
#  left_join(CM, by="DHSCLUST")
#GM_map <- unique( env_GM[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
#  left_join(GM, by="DHSCLUST")
#NG_map <- unique( env_NG[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
#  left_join(NG, by="DHSCLUST")
#ML_map <- unique( env_ML[ ,c("DHSCLUST","all_population_count_2015","aridity_2015","average_prim")]) %>%
#  left_join(ML, by="DHSCLUST")

custom_pal3 <- c(
  "1-1" = "#d3d3d3", # low x, low y
  "2-1" = "#ba8890",
  "3-1" = "#9e3547", # high x, low y
  "1-2" = "#8aa6c2",
  "2-2" = "#7a6b84", # medium x, medium y
  "3-2" = "#682a41",
  "1-3" = "#4279b0", # low x, high y
  "2-3" = "#3a4e78",
  "3-3" = "#311e3b" # high x, high y
  )


legend<-bi_legend(pal = custom_pal3,
          dim = 3,
          xlab = "Higher Climate Hazards ",
          ylab = "Higher Dropout ",
          size = 12)


bi_class_pre_h <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_pre_l <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_tem_h <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bi_class_tem_l <- c("1 - 1","2 - 1","3 - 1","1 - 2","2 - 2","3 - 2","1 - 3","2 - 3","3 - 3")
bivariate_color_scale <- data.frame(bi_class_pre_h,bi_class_pre_l,bi_class_tem_h,bi_class_tem_l,custom_pal3)


map_africa_pre_h <- env_africa %>%
    group_by(dhscc) %>%
    mutate(pre_h_class = as.numeric(cut(cum_pre_hdevi_m, 3,na.rm=T)),
           dropout_class = case_when(average_prim<0.45 ~ 3,
                                     average_prim<0.9 & average_prim>0.45 ~ 2,
                                     average_prim >=0.9 ~ 1),
           bi_class_pre_h = case_when(pre_h_class == 1 & dropout_class == 1 ~ "1 - 1",
                                      pre_h_class == 2 & dropout_class == 1 ~ "1 - 2",
                                      pre_h_class == 3 & dropout_class == 1 ~ "1 - 3",
                                      pre_h_class == 1 & dropout_class == 2 ~ "2 - 1",
                                      pre_h_class == 2 & dropout_class == 2 ~ "2 - 2",
                                      pre_h_class == 3 & dropout_class == 2 ~ "2 - 3",
                                      pre_h_class == 1 & dropout_class == 3 ~ "3 - 1",
                                      pre_h_class == 2 & dropout_class == 3 ~ "3 - 2",
                                      pre_h_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
    ungroup() %>%
    left_join(bivariate_color_scale, by = "bi_class_pre_h") %>%
  filter(all_population_count_2015>0) 
  
  
map_africa_pre_l <- env_africa %>%
  group_by(dhscc) %>%
  mutate(pre_h_class = as.numeric(cut(cum_pre_hdevi_m, 3,na.rm=T)),
         pre_l_class = as.numeric(cut(cum_pre_ldevi_m, 3,na.rm=T)),
         tem_h_class = as.numeric(cut(cum_mtemp_hdevi_m, 3,na.rm=T)),
         tem_l_class = as.numeric(cut(cum_mtemp_ldevi_m, 3,na.rm=T)),
         dropout_class = case_when(average_prim<0.45 ~ 3,
                                   average_prim<0.9 & average_prim>0.45 ~ 2,
                                   average_prim >=0.9 ~ 1),
         bi_class_pre_l = case_when(pre_l_class == 1 & dropout_class == 1 ~ "1 - 1",
                                    pre_l_class == 2 & dropout_class == 1 ~ "1 - 2",
                                    pre_l_class == 3 & dropout_class == 1 ~ "1 - 3",
                                    pre_l_class == 1 & dropout_class == 2 ~ "2 - 1",
                                    pre_l_class == 2 & dropout_class == 2 ~ "2 - 2",
                                    pre_l_class == 3 & dropout_class == 2 ~ "2 - 3",
                                    pre_l_class == 1 & dropout_class == 3 ~ "3 - 1",
                                    pre_l_class == 2 & dropout_class == 3 ~ "3 - 2",
                                    pre_l_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
  ungroup() %>%
  left_join(bivariate_color_scale, by = "bi_class_pre_l") %>%
  filter(all_population_count_2015>0) 



map_africa_tem_h <- env_africa %>%
  group_by(dhscc) %>%
  mutate(tem_h_class = as.numeric(cut(cum_mtemp_hdevi_m, 3,na.rm=T)),
         dropout_class = case_when(average_prim<0.45 ~ 3,
                                   average_prim<0.9 & average_prim>0.45 ~ 2,
                                   average_prim >=0.9 ~ 1),
         bi_class_tem_h = case_when(tem_h_class == 1 & dropout_class == 1 ~ "1 - 1",
                                    tem_h_class == 2 & dropout_class == 1 ~ "1 - 2",
                                    tem_h_class == 3 & dropout_class == 1 ~ "1 - 3",
                                    tem_h_class == 1 & dropout_class == 2 ~ "2 - 1",
                                    tem_h_class == 2 & dropout_class == 2 ~ "2 - 2",
                                    tem_h_class == 3 & dropout_class == 2 ~ "2 - 3",
                                    tem_h_class == 1 & dropout_class == 3 ~ "3 - 1",
                                    tem_h_class == 2 & dropout_class == 3 ~ "3 - 2",
                                    tem_h_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
  ungroup() %>%
  left_join(bivariate_color_scale, by = "bi_class_tem_h") %>%
  filter(all_population_count_2015>0) 


map_africa_tem_l <- env_africa %>%
  group_by(dhscc) %>%
  mutate(tem_l_class = as.numeric(cut(cum_mtemp_ldevi_m, 3,na.rm=T)),
         dropout_class = case_when(average_prim<0.45 ~ 3,
                                   average_prim<0.9 & average_prim>0.45 ~ 2,
                                   average_prim >=0.9 ~ 1),
         bi_class_tem_l = case_when(tem_l_class == 1 & dropout_class == 1 ~ "1 - 1",
                                    tem_l_class == 2 & dropout_class == 1 ~ "1 - 2",
                                    tem_l_class == 3 & dropout_class == 1 ~ "1 - 3",
                                    tem_l_class == 1 & dropout_class == 2 ~ "2 - 1",
                                    tem_l_class == 2 & dropout_class == 2 ~ "2 - 2",
                                    tem_l_class == 3 & dropout_class == 2 ~ "2 - 3",
                                    tem_l_class == 1 & dropout_class == 3 ~ "3 - 1",
                                    tem_l_class == 2 & dropout_class == 3 ~ "3 - 2",
                                    tem_l_class == 3 & dropout_class == 3 ~ "3 - 3")) %>%
  ungroup() %>%
  left_join(bivariate_color_scale, by = "bi_class_tem_l") %>%
  filter(all_population_count_2015>0)
mapview(na.omit(map_africa_pre_h), 
        zcol = "bi_class_pre_h",
        xcol = "LONGNUM", 
        ycol = "LATNUM",
        legend = F,
        crs = 4269, 
        grid = FALSE,
        #cex = 3,
        cex="all_population_count_2015",
        #alpha.regions = 0.2,
        alpha = 0.5,
        col.regions=list("#d3d3d3","#8aa6c2","#4279b0", "#ba8890","#7a6b84","#3a4e78","#9e3547","#682a41","#311e3b"))
# 
# mapview(na.omit(map_africa_pre_l), 
#         zcol = "bi_class_pre_l",
#         xcol = "LONGNUM", 
#         ycol = "LATNUM",
#         legend = F,
#         crs = 4269, 
#         grid = FALSE,
#         cex="all_population_count_2015",
#         alpha = 0.5,
#         col.regions=list("#d3d3d3","#8aa6c2","#4279b0", "#ba8890","#7a6b84","#3a4e78","#9e3547","#682a41","#311e3b"))

# mapview(na.omit(map_africa_tem_h), 
#         zcol = "bi_class_tem_h",
#         xcol = "LONGNUM", 
#         ycol = "LATNUM",
#         legend = F,
#         crs = 4269, 
#         grid = FALSE,
#         cex="all_population_count_2015",
#         alpha = 0.5,
#         col.regions=list("#d3d3d3","#8aa6c2","#4279b0", "#ba8890","#7a6b84","#3a4e78","#9e3547","#682a41","#311e3b"))
# 
# mapview(na.omit(map_africa_tem_l), 
#         zcol = "bi_class_tem_l",
#         xcol = "LONGNUM", 
#         ycol = "LATNUM",
#         legend = F,
#         crs = 4269, 
#         grid = FALSE,
#         cex="all_population_count_2015",
#         alpha = 0.5,
#         col.regions=list("#d3d3d3","#8aa6c2","#4279b0", "#ba8890","#7a6b84","#3a4e78","#9e3547","#682a41","#311e3b"))

#env_africa_single <- env_africa %>% # filter(nightlights_composite>=0) %>% # drop_na(global_human_footprint)

#mapview(env_africa_single, # zcol = “nightlights_composite”, # at = c(0,1,5,10,20,30), # xcol = “LONGNUM”, # ycol = “LATNUM”, # crs = 4269, # grid = FALSE, # cex=“all_population_count_2015”)

#mapview(env_africa_single, # zcol = “global_human_footprint”, # at = c(0,20,40,60,80,100), # xcol = “LONGNUM”, # ycol = “LATNUM”, # crs = 4269, # grid = FALSE, # cex=“all_population_count_2015”)

#mapview(na.omit(africa_map), # zcol = “aridity_2015”, # at = c(0,20,40,60), # xcol = “LONGNUM”, # ycol = “LATNUM”, # crs = 4269, # grid = FALSE, # cex=“all_population_count_2015”)

#mapview(na.omit(africa_map), # zcol = “average_prim”, # at = c(0,0.3,0.6,0.9,1), # xcol = “LONGNUM”, # ycol = “LATNUM”, # crs = 4269, # grid = FALSE, # cex=“all_population_count_2015”)

#write.csv(env_africa,“env_africa.csv”)